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We analyze the nonlinear evolution of spherically symmetric wormhole solutions coupled to a 
massless ghost scalar field using numerical methods. In a previous article we have shown that static 
wormholes with these properties are unstable with respect to linear perturbations. Here, we show 
that depending on the initial perturbation the wormholes either expand or decay to a Schwarzschild 
black hole. We estimate the time scale of the expanding solutions and the ones collapsing to a black 
hole and show that they are consistent in the regime of small perturbations with those predicted 
from perturbation theory. In the collapsing case, we also present a systematic study of the final 
black hole horizon and discuss the possibility for a luminous signal to travel from one universe to the 
other and back before the black hole forms. In the expanding case, the wormholes seem to undergo 
\ an exponential expansion, at least during the run time of our simulations. 
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I. INTRODUCTION 



Motivated by spectacular predictions such as interstellar travel and time machines [3, 0, Q there has been a lot of 
interest in wormhole physics over the last years. Although wormhole geometries are speculative because they require 
exotic matter in order to represent consistent solutions of Einstein's field equations, recent cosmological observations 
(see [ij for a review) may suggest the existence of phantom energy which could in principle be used to support a 
^Jy wormhole @, @, 0|- 

Assuming that stationary wormholes do exist within a particular matter model, a relevant question is whether or 
not they are stable with respect to perturbations. Motivated by this question, in a previous paper [8| we performed 
a stability analysis of static, spherically symmetric wormhole solutions supported by a minimally coupled massless 
ghost scalar field. As we proved, such wormhole solutions are unstable with respect to linear, spherically symmetric 
perturbations of the metric and the scalar field. More specifically, we showed that each such solution possesses a unique 
, unstable mode which grows exponentially in time, is everywhere regular and decays exponentially in the asymptotic 
region. The time scale associated to this mode is of the order of the areal radius of the wormhole throat divided by 
the speed of light. 

The aim of this paper is twofold. First, we verify the wormhole instability predicted by the linear stability analysis 
carried out in our previous paper. For this, we first construct analytically initial data which describes a nonlinear 
, perturbation of the data induced by the equilibrium configuration on a static slice. The time evolution of this data 
is obtained by numerically integrating the nonlinear field equations and we show that small initial perturbations do 
• i-h , not decay to zero but grow in time and eventually kick the wormhole throat away from its equilibrium configuration. 
The growth of the departure of the areal radius of the throat from its equilibrium value as a function of proper time 
at the throat is computed for different parameters of the initial perturbation and shown to agree well with the time 
■ - - ■ scale predicted by the perturbation calculation. 

Our second goal is to determine the end state of the nonlinear evolution. We find that depending on the details 
of the initial perturbation the wormhole either starts expanding or collapsing. In the collapsing case, we observe the 
formation of an apparent horizon. By computing geometric quantities at the apparent horizon, we obtain a strong 
indication that the apparent horizon eventually settles down to the event horizon of a Schwarzschild black hole. By 
analyzing the convergence of light rays near the apparent horizon we obtain an estimate for the location of the event 
horizon and find that although asymptotically it seems to agree with the location of the apparent horizon it lies inside 
it at early times. Despite the rapid collapse, we find that it is possible for a photon to travel from one universe to the 
other and back without falling into the black hole. In the expanding case, we do not observe any apparent horizons 
and our simulations suggest that the wormhole expands forever. In particular, for wormholes which are reflection 
symmetric we find that the areal radius at the point of reflection symmetry grows exponentially as a function of proper 
time. Furthermore, we point out some problems with defining the wormhole throat in an invariant way. Specifically, 
we show that the definition of the throat as a minimum of the areal radius over each time-slice does not yield an 
invariant three surface, and we give an example showing that this three surface depends on the slicing condition. 

Furthermore, we perform a careful study for the relation between the radius of the apparent horizon of the final 
black hole and the amplitude of the perturbation. By varying this amplitude, three regions are found: i) positive 
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values of the amplitude result in the collapse to a black hole, ii) negative values of the amplitude lying sufficiently 
close to zero trigger the explosion of the solution, iii) a second threshold for the amplitude is found on the negative 
real axis below which the wormholes collapse to a black hole. 

This paper is organized as follows. In section [II] we formulate the Cauchy problem for obtaining the nonlinear time 
evolution of spherically symmetric configurations, review the static wormhole solutions and describe our method for 
obtaining initial perturbations thereof. Our numerical implementation and the diagnostics we use in order to test 
our code and analyze the data obtained from the evolution are described in section IIII1 In section IIVI we analyze the 
collapsing case and indicate that the final state is a Schwarzschild black hole. By analyzing the resulting spacetime 
diagram we also mention different scenarios for a light ray to travel from one universe to the other and back. The 
expanding case is analyzed in section [V] and the study for the relation between the final state and the amplitude of 
the perturbation is carried out in section fVTl Finally, we summarize our findings and draw conclusions in section fVIII 

Some of the conclusions presented in this paper have been presented in Ref. Q for the wormholes with zero ADM 
masses based on a numerical method using null coordinates. 



II. THE CAUCHY PROBLEM FOR TIME-DEPENDENT SOLUTIONS 



In this section we formulate Einstein's equations coupled to a massless ghost scalar field in a suitable form for the 
numerical integration of spherically symmetric wormhole solutions. We use geometrized units for which the speed of 
light and Newton's constant are equal to one. We write the metric in the form 

ds 2 = -e 2d dt 2 + e 2a dx 2 + e 2c (dd 2 + sin 2 ■& dip 2 ) , (1) 

where the functions d = d(t,x), a = a{t,x) and c = c(t,x) depend only on the time coordinate t and the radial 
coordinate x. Similarly, we assume that the scalar field $ = <E>(£, x) only depends on t and x. The coordinate x ranges 
over the whole real line, where the two regions x — > oo and x — ► — oo describe the two asymptotically flat ends. In 
particular, we require that the two-manifold (M,g) — (R 2 ,—e 2d dt 2 + e 2a dx 2 ) is regular and asymptotically flat at 
x — > ±00, and that the areal radius r = e c is strictly positive and proportional to |x| for large For this reason, it 
is convenient to replace c by c where 

c = c+Uog{x 2 + b 2 ) (2) 

and b > is a parameter. Next, we find it convenient to choose the following family of gauge conditions: 

d = a + 2Ac, (3) 

which relates the lapse e d to the functions a and c parametrizing the three-metric. Here, A is a parameter that may 
take any real value in principle, but for our simulations we will consider the choices A = and A = 1. The metric now 
takes the form 

ds 2 = e 2a (_ e 4Ac^2 + ^2) + ^2 + tf^Tc ^2 + ^2 # ^2) _ (4) 



A. Evolution and constraint equations 

Einstein's equations yield the Hamiltonian constraint Ti := —e a ~ d Gu + K[e a ~ d & 2 -I- e d ~ a Q 2 ]/2 — 0, the momentum 
constraint M := — [R x t~ ^t^x) — and the evolution equations R$$ = R vv = and e d ~ a (R xx — k$ 2 ) + (1+\)H = 0, 
where here and R^ refer, respectively, to the components of the Einstein and Ricci tensor, k = —8ir and $ t := <9 t $ 
and (fr x :— d x &. These evolution equations, together with the wave equation for $ yield the following coupled system 
of nonlinear wave equations for the quantities a, c and $: 



d t (e- 2X£ a t ) - d x 



,2Ac 



2A£ 



b 2 



e -2\c ^ + A ) g2 + 2 \a t c t ] + e 2Ae [(1 + 3A)c 2 - 2Xc x (a x + 2Xc x )] 



2a+2(A-l)c 

-(1 + A) - ; , 9 + J [(A + l) e - 2Ae $ 2 + (A - l)e 2A5 $ 2 ] = 0, 



dt (e 2(1 - 



b 2 
1 



b 2 



d x 



b 2 



e 2 ^ s ((x 2 + b 2 )c x +x) 



(x 2 + b 2 )e 2 ^ 1+x ^ a 



e 2{a+Xc) 

x 2 + b 2 



= 0, 



(5) 

(6) 
(7) 
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which is subject to the constraints 

2a+2(A-l)c 

U = e 2X5 (2c xx + 3c 2 x - 2a x c x ) - e - 2Xs c t (c t + 2a t ) - \ 2 + ^ + ~ [e~ 2X ^ 2 + e 2XE <S> 2 x ] = 0, (8) 

M = 2c tx + 2c t (c x - 2Xc x - a x ) - 2a t c x + k $> t $ x = 0. (9) 

For A = we see that the principal part of the evolution equations is given by the flat wave operator df — d 2 . 
Therefore, the characteristic lines coincide with the radial null rays which for A = are given by the straight lines 
t ± x = const. Therefore, a property of the gauge condition ([3]) with A = is that it cannot lead to shock formations 
due to the crossing of characteristics. We found this gauge to be convenient for studying expanding wormholes. For 
the collapsing case, on the other hand, we found the gauge © with A = 1 more useful since it seems to be avoiding 
the singularity after the black hole forms. 



B. Propagation of the constraints 



It can be shown that the Bianchi identities and the evolution equations (|5l6l7p imply that the constraint variables 
Ti and M. satisfy a linear evolution system of the form 

d t H = e 2X5 d x M+l.o., (10) 
d t M = e 2X5 d x H + l.o., (11) 

where l.o. are lower order terms which depend only on Ti. and M. but not their derivatives. For a solution on the 
unbounded domain — oo < x < oo with initial data satisfying 7i = M, = it follows that the constraints are 
automatically satisfied everywhere and for each t > 0. Therefore, it is sufficient to solve the constraints initially. 
If artificial time-like boundaries are imposed, this statement is still true provided suitable boundary conditions are 
specified. For example, imposing the momentum constraint M — at the boundary ensures that the constraints 
remain satisfied if so initially. 

C. Boundary conditions 

Since we have three wave equations, at the artificial boundaries we apply outgoing wave boundary conditions for 
the three fields a, c and We assume that these fields behave like spherical waves far away from the origin, i.e. if 
x > 0: 

/(*,*) = ^^, (12) 
x 

where v — e 2Xc is the speed of propagation. In order to apply the boundary conditions, we use the following equation, 

-d t f + d x f+£ = 0. (13) 
v x 

We use a similar procedure when x < 0. It is clear that such conditions do not preserve the constraints and 
consequently, the solution will be contaminated with a constraint-violating pulse traveling inwards in the numerical 
domain. In order to avoid the contamination of the analyzed data, we push the boundaries far from the throat such 
that the region where we extract physics is causally disconnected from the boundaries. 

D. Static wormhole solutions 

In the static case, all wormhole solutions can be found analytically [l(| [HI G3 and are given by the following 
expressions 

d = — a = 71 arctan ^— ^ , (14) 

c = —71 arctan ^— J , (15) 

$ = $1 arctan (J^J , (16) 
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where the parameters $1 and 71 are subject to the condition = 2(1 + jf). The wormhole throat is located at 

Xthroat = lib and has areal radius rthroat = &\A + 7ie~ 71 arctan ('w). The ADM masses at the two asymptotically flat 
ends x — + ±00 are = 671 exp(— jin/2) and m_oo = — 671 exp(7i7r/2), respectively (see Ref. Q for a derivation 
and details). In particular, 71 = yields wormhole solutions with zero ADM mass at both ends while in all other 
cases the ADM masses at the two ends have opposite signs. 



E. Initial data describing a perturbed wormhole 



In order to study the nonlinear stability of the static solutions described by equations (|14I15I16[) we construct initial 
data corresponding to an initial perturbation of a static solution. For simplicity, we assume that the initial slice is 
time-symmetric in which case the momentum constraint M. = is automatically satisfied. Therefore, we only need 
to solve the Hamiltonian constraint which, for a t — ct = = simplifies to 

2x 1 — e 2 ( a ~ e ) b 2 k 

2c xx + 3c* - 2a x c x + -^{^ ~ «.) + ^ + = ~ ^ ^ 

Since the unperturbed wormhole solution is known in analytic form a simple way of obtaining initial data representing 
a perturbation thereof is to perturb the metric quantities a and c by hand in such a way that the left-hand side of 
equation (TTT)) is nonnegative and to solve equation (|17p for Q x . Since the left-hand side of Equation (fTT|) is everywhere 
positive for the analytic solution, it remains positive at least for small enough perturbations. Moreover, it is clear 
that any time-symmetric and spherically symmetric initial data can be obtained by this method. Next, we notice 
that it is sufficient to consider the case where the metric component a is unperturbed. Indeed, a perturbation of a 
may be absorbed by a redefinition of the coordinate x which, at the physical level, does not change the initial data 
and its Cauchy development. Therefore, it is sufficient to perturb the quantity c which is related to the areal radius. 
For this work, we choose to perturb it with a Gaussian pulse. More precisely, we choose 



^pcrt — ^static ) 

(18) 

Cpcrt — ^static ~t~ ^ ^ ^ c j (-^) 



where e c is the amplitude, x c the center and <r c is related to the width w at half maximum of the pulse through 
w = 2 -y/log 2 a c . As will be shown later, different values of e c and a c produce two different scenarios: the collapse of 
the wormhole to a black hole or a rapid expansion. Notice that because of the exponential decay of the perturbation 
as |x| — > ±00, the ADM masses of the perturbed solution is equal to the ADM masses of the unperturbed, static 
solution. 



III. IMPLEMENTATION AND DIAGNOSTICS 



In this section we describe our numerical method and different tools used to analyze the data obtained from 
the simulations. The numerical method is based on a second order centered finite differences approximation of the 
evolution equations ([S][7]) and the constraint equations (|8I9[) . We only perform the evolution on a finite domain with 
artificial boundaries at a finite value of x where we implement the radiative-type boundary condition (|13p with an 
accuracy of second order for the evolution variables a, c and <I>. A method of lines using the third order Runge-Kutta 
integrator is adopted. Throughout the evolution, we monitor the constraint variables 7i and M. and check that they 
converge to zero as resolution is increased. 

From now on we rescale the coordinate x in such a way that 6=1. 



A. Unperturbed wormhole 



As a test, using the implementation aforementioned we perform a series of simulations for the unperturbed, static 
wormholes. We start a simulation with the massless case 71 = 0, for which it is expected that the fields a, c, $ 
remain time- independent. However, discretization errors are enough to trigger a non-trivial time dependence of these 
functions. We have verified that when increasing the resolution, the numerical error remains small for a larger period 
of time and the numerical solution converges to the exact, static solution with second order. In figure [1] we show the 
convergence of the constraint variable TL for a wormhole using the gauge parameter A = 1. 



5 




18-07 1 ' 1 ' 1 

5 10 15 20 

t 



FIG. 1: The L2 norm of the Hamiltonian constraint variable TL as a function of coordinate time for an unperturbed wormhole 
with gauge parameter A = 1. Different resolutions are shown. As is apparent from the plot at each fixed time, the error 
decreases with increasing resolution. 



B. Extremal and marginally trapped surfaces 

Let us consider at — const hypersurface £4. There are two different types of interesting two-surfaces in this 
hypersurface: The first are extremal surfaces which consist of two-surfaces with a local extrema of their area. If s l 
denotes the unit normal to this surface, this means that the divergence DiS 1 vanishes, where D refers to the covariant 
derivative associated to the induced three- metric on E t . The second type of two-surfaces of interest are marginally 
trapped surfaces. Denoting by l a and k a the future-pointing outward and inward null-vectors orthogonal to this 
two-surface, a marginally trapped surface is defined by the requirement that the expansion 9- along k a is strictly 
negative while the expansion 9+ along l a vanishes. In terms of the induced metric 7 Q/ 3 on the two-surface, this means 
that 

0- := 7a/ 3 V a fe /3 < 0, 9+ := lal) V a lP = 0. 

Notice that a rescaling of k a or l a by a positive function does not affect this definition. On the other hand, for 
wormhole topologies the notions "outward" and "inward" are observer-dependent: they depend on which asymptotic 
end these vectors are viewed from. 

With respect to t = const slices of the spherically symmetric metric (JTJ we have s t di = e~ a d x> k a d a = e~ d d t — e~ a d x , 
l a d a = e~ d dt + e~ a d x (with respect to the asymptotic end x — ► 00), and 

D lS l = 2e~ a c x , 0_ = 2(e- d c t - e~ a c x ), 9+ = 2(e- d c t + e~ a c x ). (20) 

Therefore, an extremal surface is a sphere S 2 for which c x = 0, and a marginally trapped surface a sphere for which 
e~ d Ct = —e~ a c x and c x > 0. In particular, according to our definition, an extremal surface cannot be a marginally 
trapped surface. Finally, we define a throat to be an extremal surface which (out of the many extremal surfaces 
that might exist) is one with least area. In a similar way, an apparent horizon is an outermost marginally trapped 
surface. Notice that these concepts are not defined in a geometrically invariant way. For example, it is known [l3l | that 
the Schwarzschild spacetime possesses Cauchy surfaces which do not contain any trapped surfaces and nevertheless 
come arbitrarily close to the singularity. Similarly, we will provide a numerical example below that shows that the 
three-surface obtained by piling up the throats in each time slice is not invariant but depends on the time foliation. 

C. Geometric quantities 

Throughout evolution, we monitor the values of the following geometric quantities: 

r = e c = ^/x 2 + b 2 e e , (21) 
M = r -(l - gf»V ar ■ V b r) = | [l - e 2 %~e' 2d c 2 + e " a «<£)] , (22) 
L = <f b V Q $ • V fc $ = -e- 2d <5> 2 t +e- 2a <f> 2 x , (23) 
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where quantities with a tilde refer to the two-metric g = —e 2d dt 2 + e 2a dx 2 . As mentioned above, r — r(t,x) is the 
areal radius of the sphere at constant t and x. M is the Misner-Sharp mass function [l4[ and L is the norm of the 
gradient of the scalar field. In particular, the knowledge of these functions allows the computation of the following 
curvature scalars: The Ricci scalar R associated to the two-metric g, the Kretschmann scalar associated to the full 
metric I := R al3lS R af3l s, and the square of the Ricci tensor, J := WR^. Using Einstein's field equations one 
obtains the following expressions, 

4M 

R = ^- + kL, (24) 



A&M 2 8kM T „ , , 
1 = ~^ + ^T L + 2k L > ( 25 ) 



r 

J = k 2 L 2 . (26) 



D. The construction of conformal coordinates 

The gauge choice © with A = has the property of yielding conformaUy flat coordinates (t, x) for the two-metric 
g = —e 2d dt 2 + e 2a dx 2 . In these coordinates, the null rays are simply given by the straight lines t ± x = const. 

On the other hand, in many simulations, the choice A = 1 seems to work better than the choice A = 0. In order 
to compare the coordinates (T,X), say, constructed with the gauge choice A = 1 to the conformally flat coordinates 
(i, x) obtained with A = one can proceed as follows. Suppose we are given to us a two-metric g = g a bdx a dx b with 
signature (— 1, 1). We are interested in finding local coordinates (t, x) such that 

g = e 2a {-dt 2 + dx 2 ), (27) 

with a conformal factor e a . We first note that if i a b denotes the volume element associated to g, the coordinates t 
and x must satisfy 

V a t = ±e a b V b x , (28) 

since by equation (|27p V ' a t and V x are orthogonal to each other and their norms are equal in magnitude. It follows 
from equation {28| that t and x must also satisfy the wave equation 

V Q V a i = V a V Q ir = 0. (29) 

Therefore, in order to construct the coordinates (t,x), we have to solve equations (|2T)|) subject to the constraint 
(|28p. The resulting functions (t,x) indeed satisfy (j2"7| since 

-9 U = -g ab (V a t)(y b t) = .g ah (V a x)(V,x) = ~ g * x , 

gte = gF* at )(y b x)=0. 

The conformal factor is then obtained from e a — 1/ ' \Jg xx '. Notice that equation (|29|) implies that both C a := 
X7 a t ± e a b VbX and D a := e a b Cb are conserved, i.e. V a C = X7 a D a = 0. Therefore, it is sufficient to solve the 
constraint C a = on a Cauchy slice. 

For example, suppose that g — e 2A [~e 4c dT 2 + dX 2 } as is the case for our simulations with A = 1. Then, equa- 



tions ([29]) yield 

d T (e- 2s d T t) = d x (e 2£ d x t), 
d T (e' 2s d T x) = d x (e 25 d x x), 
and the initial data for the functions t = t(T, X) and x = x(T, X) satisfies 

t(0,X) = 0, 

x(0,X) = X, 

d T t(0,X) = e 2E d x x(0,X) = e 25(0 ' X \ 

d T x(0,X) = e 2£ d x t{0,X) = 0. 



The last two equations follow from equation (f28|) where we have chosen the + sign such that t and T point in the 
same direction. The so obtained functions t = t(T, X) and x = x(T, X) can be used, for instance, to plot a T = const 
surface in the t — x diagram. They may also be useful to find the null rays which are given by the lines with constant 
t + x or t — x. 
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IV. COLLAPSE TO A BLACK HOLE 

For simplicity, we start with the evolution of a perturbed, massless wormhole with the perturbation consisting of 
a Gaussian pulse as in (fH)| with x c = so that the perturbation is centered at the throat. We start with positive 
values of the amplitude e c and find that such perturbations induce a collapse of the throat. The typical behaviour of 
such a collapsing wormhole is presented in figures [2] to [9j 

In figure[5]we show the convergence of the Hamiltonian constraint as a function of coordinate time and the evolution 
of the areal radius which exhibits the collapse of the wormhole throat. 




FIG. 2: (Left panel) L2 norm of the Hamiltonian constraint variable H as a function of coordinate time for a perturbed 
wormhole with parameters e c — 0.001 and a c = 0.5 and gauge parameter A = 1. Various resolutions are shown. We have 
checked that these constraint variables converge to zero with second order. (Right panel) Areal radius as function of the 
coordinate x at different times. In this scenario the value of the throat's areal radius starts at rthroat = 1 and collapses to zero 
in finite coordinate and proper time (see also figure [3] below). 



A. Time scale of the collapse 

In order to check the consistency of our numerical calculations with the time scale associated to the unstable mode 
found in our linear stability analysis, we perform a systematic study of the speed of the collapse of the throat. For 
simplicity, we only consider massless wormholes with a centered perturbation. Our procedure for quantifying the 
speed of the collapse is the following: 

i) Compute at each time the throat's location by finding the minimum of the areal radius. For the reflection 
symmetric case, this minimum is located at x = 0. 

ii) Compute the proper time r integrating the lapse function e d at the throat, 
hi) Plot the throat's areal radius as a function of proper time. 

iv) Fit the throat's areal radius as a function of proper time to the function r(r) = 1 — e^^ T ^ Pl ^ P2 , from the initial 
radius of the throat until a chosen minimum value of r, called r cut . We interpret the parameter P2 as the time 
scale of the solution. The reason to cut off the values of r is that we want to compare P2 with the results from 
perturbation theory which are expected to be valid as long as the departure from the equilibrium configuration 
is small. 

v) Repeat the above analysis for fixed initial parameters x c — and a c — 0.5 and several values of the amplitude 

In figure [3] we show the results for the time scale and its dependency on the cut-off value r cut , the initial parameter 
e c and the resolution used for the numerical evolutions. Two main features are: 1) The time scale for small values of 
e c and r cut < 1 approach the one calculated from linear perturbation theory Q, 2) for large values of e c and small 
values of r cut , the time scale is always smaller than the result of perturbation theory, indicating that nonlinear terms 
tend to accelerate the collapse. 
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FIG. 3: (a) We show a typical evolution of the throat's areal radius r t ; lroat in terms of proper time r for the case of a collapse 
induced by a centered perturbation. The inset shows the same function with a logarithmic scale for r t hroat and illustrates how 
rthroat converges to zero in finite proper time r. Notice also that r freezes with respect to coordinate time which indicates 
that the lapse collapses to zero near x — 0. (b) The time scale p2 is calculated for different values of e c and various resolutions 
using r cu t = 0.9. Second order convergence of the time scale for each value of the amplitude is manifest. We predict the value 
at infinite resolution using Richardson extrapolation. For this particular case, the time scale for the smallest perturbation 
(e c = 0.001) shows the highest time scale of all cases, namely 0.836. This value is comparable to the one obtained using linear 
perturbation theory which is 0.846, see Table I in Ref. [|. (c) to (f) The same analysis for different values of r cu t. These 
plots indicate a decreasing monotonic behaviour of the time scale as r cu t decreases, indicating that nonlinear terms tend to 
accelerate the collapse. 
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B. Formation of an apparent horizon 

During the collapse, we observe the formation of an apparent horizon. This is shown in the left panel of figure |4l 
where we also show a bundle of outgoing null rays whose intersection points with the initial slice are fine-tuned in 
such a way that these rays lie as close as possible to the apparent horizon for late times. As the plot suggests, this 
bundle of light rays separates photons that might escape to future null infinity from photons that are trapped inside 
a region of small x. Therefore, we expect this bundle to represent a good approximation for the event horizon of a 
black hole. As a consequence, our simulations give a strong indication for the collapse of the wormhole into a black 
hole. 

The bold black line in the left panel of figure [4] indicates a possible trajectory of a photon traveling from one universe 
to the other and back, allowing it to reach the asymptotic region of the home universe without getting trapped by 
the black hole. In the right panel the future null directions at the apparent horizon location are presented, showing 
that the latter is a time-like surface. 




FIG. 4: These plots correspond to the previously shown collapsing case. (Left panel) We show the apparent horizon surface 
and a bundle of null rays which depart from a surface identified with an approximation of an event horizon. The portion of 
the apparent horizon for x > (x < 0) corresponds to the outermost marginally trapped surface as viewed from an observer 
at x — ► oo (x — » — oo). Correspondingly, the portion of the bundle of null rays emanating near x = —6 (x = +6) represents the 
event horizon with respect to the universe x > (x < 0). With respect to an observer in either universe, it can be seen that 
the apparent horizon lies outside the event horizon. The bold line indicates a possible trajectory of a null geodesic going from 
one universe to the other and reflected back. (Right panel) We show the future null directions at the location of the apparent 
horizon for x > 0; the arrows departing from the apparent horizon indicate the in- and outgoing directions of light cones. This 
result indicates that the apparent horizon is actually a time-like surface. 

Somehow unusual features of the collapse are that the apparent horizon is not contained in the black hole region 
and that the areal radius of both the apparent and event horizons decreases in time. This is illustrated in figure OH 
where we plot the trajectory of the apparent horizon and the bundle of outgoing null rays in the t — r diagram. Notice 
that for a spacetime satisfying the null energy condition and cosmic censorship, propositions 9.2.1 in 1 1 51 ] and 12.2.3 
in [l6| show that an apparent horizon is contained in the black hole region. Under the same assumptions theorem 
12.2.6 in [ll| and [U [H[ establish that the area of the event horizon cannot decrease in time. In the present case, 
these results do not apply since the null energy condition is not satisfied. The decrease of the area of the horizons 
could be related to the absorption of the ghost scalar field by the black hole. 

C. The final state 

In order to analyze the final state of the collapsing wormholes we compute the scalars L and K = ^f- from 
equations (|22|23|) at the apparent horizon location. At the event horizon of a Schwarzschild black hole these scalars 
are zero and one, respectively. As shown in figure El these quantities do indeed converge to the corresponding 
Schwarzschild values for large times and high spatial resolutions, indicating that the final state is a Schwarzschild 
black hole, at least in the vicinity of the apparent horizon. 

We try different initial parameters and perform runs with centered and off-centered perturbations, and in all cases 
where the wormhole starts collapsing the result is the same: an apparent horizon forms, and it seems to settle down 
to the event horizon of a Schwarzschild black hole. We also investigate the mass of the final black hole as a function of 
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FIG. 5: The trajectory of the apparent horizon and the bundle of outgoing null rays in the t — r diagram. Notice that the areal 
radius of both the apparent and the event horizons decreases in time. Even though these trajectories cross at t ~ 5 the event 
horizon lies inside the apparent horizon at all times, as is apparent from the t — x diagram shown in figure [4] above. 
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FIG. 6: The scalars L and K evaluated at the apparent horizon versus coordinate time for different resolutions. For high 
resolutions and late times, Lah converges to zero and Kah to one. This indicates that the apparent horizon converges to 
the event horizon of a Schwarzschild black hole. The collapse is triggered by the off-centered perturbation with e c = 0.00125, 
a c = 0.5, x c = 2. 

the parameters of the perturbation. In figure [7] we present the evolution of the apparent horizon radius as a function 
of coordinate time for small, positive values for the amplitude e c and different values for the width a c . For such values 
it seems that the final black hole mass is universal and given by ttiah ~ 0.22. The dependency of the black hole's 
final areal radius for a wider range of initial parameters will be analyzed in section IVT1 

D. The distribution of the scalar field 

In order to investigate the time evolution of the distribution of the scalar field, we show in figure [5] snapshots at 
different times of L along the spatial domain for the collapse of a zero mass configuration. A pulse of scalar field 
departs from the location of the apparent horizon region and propagates outwards with negative values. Since we 
know that the ADM mass of the spacetime is zero and that the mass of the apparent horizon is positive, we expect 
these pulses of scalar field to radiate negative energy. 

Finally, we also analyzed the collapse of massive wormholes which are characterized by a positive value of the 
parameter 71. The results are qualitatively similar to those found for the zero mass case. For completeness, we show 
the behaviour of the areal radius for a collapsing case and the scalars calculated at the apparent horizon in figure [HI 
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FIG. 7: Apparent horizon radius as function of the perturbations. (Left) The width of the perturbation is fixed to a c — 0.5 
and the amplitude of the perturbation ranges from e c = 0.000 to e c = 0.010 in intervals of Ae c = 0.001. (Right) The amplitude 
of the perturbation is fixed to the value e c = 0.002 and the width varies from a c — 0.1 to cr c = 1.0 in intervals of Aa c = 0.1. 
These plots suggest that the final value for the areal radius of the apparent horizon is independent of the perturbation in this 
range of parameter space. 
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FIG. 8: Snapshots of L for a zero mass collapsing configuration and a centered perturbation. Notice the outgoing pulse with 
negative amplitude developing outside the apparent horizon region. 



V. THE EXPANDING CASE 



Perturbations that do not result in a collapse to a black hole induce a rapid growth of the wormhole throat. In 
order to analyze this case, we performed a series of simulations for the perturbed wormhole solutions starting with 
the massless case 71 = 0. In figure [TOl we consider the evolution of a reflection symmetric perturbation and show the 
areal radius of the spheres at the points of reflection symmetry x — as a function of proper time. The exponential 
growth of this function indicates a rapid growth of the wormhole. 
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FIG. 9: (Left panel) Snapshots of the areal radius in terms of x for several times. The particular behaviour of the massive case 
is characterized by an initial boost of the solution in the coordinates we use. This simulation was carried out using the physical 
parameter 71 = 0.015 and the perturbation parameters e c = 0.002, a c — 0.5 and x c — 0. (Right panel) The values of K and L 
evaluated at the apparent horizon. Once again we find that these scalars converge to their Schwarzschild values for late times 
and high resolutions. 
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FIG. 10: (Left panel) We show a typical evolution of the throat's areal radius versus proper time t for the expanding case 
under the action of a centered perturbation. (Right panel) The time scale P2 is calculated for different values of e c and various 
resolutions using r cu t = 1.1. Second order convergence of the time scale for each value of the amplitude is manifest and the 
Richardson extrapolation value is shown. The time scale for the smallest perturbation lies near 0.85. 



A. Time scale of the expansion 

Also shown in figure [10] is the time scale P2 associated to the exponential growth for different resolutions and 
different values of the initial amplitude e c of the perturbation. The function used to fit the radius of the throat in this 
case is r(r) = 1 + e( T ~ Pl ^ P2 . In the limit e c — > the time scale is estimated to be 0.85 which is in good agreement 
with the prediction from perturbation theory, see Table I in Ref. [1]. 

B. Slice-dependence of the throat 

Next, we performed runs with the two gauge values A = and A = 1 with perturbations centered at the throat. 
While we did not find the formation of any apparent horizons in any of these two gauges, there is an interesting 
difference in the behaviour of the areal radius r as a function of the x coordinate. This is shown in figure 111! For 
A = this function always has its minimum at x = 0, and one is tempted to define the throat's location to be at 
x = 0. However, looking at the results for A = 1, one finds that although initially the minimum of r as a function of 
x lies at x — 0, at late times the function r possesses a maximum at x — and two minima, one at a negative value 
of x and the other at a positive value of x. Therefore, from the results with A = 1, one is led to the conclusion that 
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two throats develop. 
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FIG. 11: Areal radius r versus x for A = and A = 1 for the case of expansion using a centered perturbation. When A = 0, the 
areal radius grows monotonically and preserves the notion of a single throat in coordinate time. In the case A = 1 the areal 
radius develops a sort of two-throat image as can be seen in the inset region of the plot. 

It turns out that this apparent paradox is an effect of the different slicing conditions resulting from the gauge 
condition ([3]) with A = and A = 1, respectively. In order to explain this, denote by (i, x) the conformal coordinates 
corresponding to the evolution with A = 0, and let (T, X) denote the coordinates obtained in the evolution with A = 1. 
The relation between these two coordinate systems can be obtained from the method described in section IIII Dl In 
figure[12]we show a spacetime diagram based on the conformal coordinates (t, x) where we plot the surfaces of constant 
areal radius r and the space-like slices of constant T. As can be seen, at late enough times, the surfaces of constant 
T may cross the lines of constant areal radius four times. Therefore, as one moves along a T — const slice from 
one asymptotic end towards the other, the areal radius decreases until a minimum is reached, then the areal radius 
increases again until a local maximum is reached at x — 0; then r decreases again until a local minimum is reached, 
and finally the areal radius increases again. Therefore, the two-throat image is a pure gauge effect, and the definition 
of the throat's location as a local minimum of the areal radius over a given time slice is not defined in a geometrical 
invariant way, but strongly depends on the time foliation. 




FIG. 12: Surfaces of constant T (the time coordinate in the coordinate system obtained from the evolution with A = 1) and 
constant areal radius r. As can be seen, a T = const surface may cross a r — const surface four times, indicating that the 
two-throat image is a pure gauge effect. 



C. The distribution of the scalar field 

In order to investigate the time evolution of the distribution of the scalar field for the expanding case we show 
in figure [TJJ] snapshots at different times t of the scalar L along the spatial domain for the expansion of a zero mass 
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wormhole. As can be seen from this plot, the amplitude of this quantity increases in time while its shape remains 
essentially the same. Therefore, in contrast to the collapsing case, the scalar field does not seem to dissipate away 
but actually gains in strength near the point of reflection symmetry x — 0. 
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FIG. 13: Snapshots of L for different times t for the case of an expanding wormhole with zero mass and centered perturbation. 
Finally, we also perform an evolution for a massive wormhole with 71 = 0.015. This is shown in figure [T4l 





FIG. 14: Areal radius r versus x for A = and A = 1 for a massive case. As in the collapsing case, the throat shifts along the 
x coordinate and then the areal radius expands. The difference for different values of A is also manifest. 



VI. DEPENDENCY OF THE FINAL STATE ON THE INITIAL PERTURBATION 

We want to analyze the dependency of the final state of the system on the initial perturbation. As we mentioned 
before, the two possible final states of the evolution are the collapse into a black hole or the expansion of the throat. 
In the former case, we can measure the radius of the apparent horizon of the final black hole. In the latter, there 
is no formation of black holes. We perform a systematic analysis on the dependency of the apparent horizon's areal 
radius as a function of the initial perturbation. As before, we perturb the function c with a Gaussian pulse centered 
at xq — with fixed width a c — 0.5, varying the amplitude e c from e c ~ +0.1 to e c ~ —0.1 in intervals of Ae c = 0.01. 



15 



When the values of the amplitude are positive the system collapses into a black hole. We found that the radius of 
the black hole decreases from tah = 0.51 for e c = 0.07 to tah — 0.45 for e c = 0. Once the values of the amplitude 
become negative the behaviour of the evolution changes. In the interval e c — [—0.04, 0.00) the wormhole expands 
instead of collapsing and no apparent horizon is found. It is tempting to stop here and assume that this behaviour 
will continue. Nevertheless, if we keep decreasing the value of the amplitude we find that for values in the interval 
e c = [—0.09, —0.05] the system collapses again into a black hole. In this range, the radius of the black hole increases 
from tah = 0.49 to Tah = 0.53. These results are presented in figure [TBI for three different resolutions in order to 
emphasize the convergence of the apparent horizon's areal radius. 
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FIG. 15: The radius of the apparent horizon is shown for various values of the amplitude of the perturbation. There are three 
regions in the plot: i) one on the left, where a finite mass horizon forms, ii) a second region in the middle for which there is 
expansion and iii) another region for positive e c for which there is collapse again. For completeness, we also show these results 
for various resolutions which indicate nearly second order convergence at each point of the plot. The four points with tah = 
correspond to expanding wormholes in which case there is no apparent horizon. 



VII. CONCLUSIONS 



Based on numerical methods, we analyzed in this work the nonlinear stability of static, spherically symmetric 
general relativistic wormhole solutions sourced by a massless ghost scalar field. This complements the linear stability 
analysis performed in Q where we proved that all such solutions are unstable with respect to linear fluctuations. 
In particular, our numerical simulations confirm the instability predicted by linear theory and show that static, 
spherically symmetric wormholes sourced by a massless ghost scalar field are also unstable with respect to nonlinear 
fluctuations. Furthermore, we have checked that the time scale associated to the instability agrees with the one 
computed from perturbation theory, at least when computed over times for which the departure from the equilibrium 
configuration is small. 

Our numerical simulations also reveal that depending on the initial perturbation, the wormholes either collapse 
to a Schwarzschild black hole or undergo a rapid expansion. This confirms the results in in the zero mass case 
and shows that a similar result holds for massive wormholes. In the collapsing case, we reach this conclusion by first 
observing the formation of an apparent horizon whose areal radius converges to a fixed positive value at late times. 
By computing geometrical quantities at the apparent horizon and analyzing the behaviour of null geodesies in the 
vicinity of the apparent horizon at late times we are led to the conclusion that the final state of the collapse is a 
Schwarzschild black hole. Further clues in support of this scenario are provided by the observation that the scalar 
field disperses. 

In the expanding case, the wormhole starts growing rapidly. For reflection symmetric wormholes we find that the 
areal radius of the spheres at the points of reflection symmetry grow exponentially with respect to proper time. On 
the other hand, we have also found that the definition of the throat as the three-surface obtained by determining the 
global minimum of the areal radius in each time slice is problematic in the sense that it may depend on the time 
foliation. 

We have also performed a systematic study of the dependency of the final state on the initial parameters of 
the perturbation. Using only reflection symmetry initial data with a Gaussian perturbation profile we investigated 
the possible universality of the final black hole's mass. We found that for small, positive amplitudes of the initial 
perturbation the final mass does not vary much with respect to the initial width of the perturbation. For initial 
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amplitudes which are negative and small in magnitude, on the other hand, the wormholcs expand, the limiting case 
of zero amplitudes corresponding to a threshold in parameter space separating collapsing wormholcs from expanding 
ones. By exploring a wide region of negative values for the initial amplitude we also found a second threshold below 
which the wormholes collapse again into a black hole. 

Since a tiny perturbation may cause the wormhole to collapse into a black hole, it is unlikely that the wormholcs 
we have analyzed in this article are useful for interstellar travel or building time machines. Furthermore, the time 
scale associated to the collapse is of the order of the throat's areal radius divided by the speed of light, which is of the 
order of a few microseconds for a throat with areal radius of one kilometer. We have analyzed the possible scenarios 
in which a test particle may use the wormhole to tunnel from one universe to the other and back before the black 
hole forms. As a consequence of the instability, an observer moving on the path of such a test particle cannot explore 
an arbitrarily large region of the other universe if he wants to travel back home. 
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